#include <stdio.h>    // 标准输入输出库，用于打印结果
#include <math.h>     // 数学函数库，用于ceil和sqrt计算
#include <cuda.h>     // CUDA核心库，支持GPU编程

typedef float FLOAT;  // 定义浮点类型别名，方便后续切换为double
#define USE_UNIX 1    // 使用Unix时间函数（1为启用，0为禁用Windows时间函数）

/* 获取块ID：2D网格 */
// 将二维块索引转换为一维序号，公式：块ID = x轴块索引 + y轴块索引 * x轴块数量
#define get_bid() (blockIdx.x + blockIdx.y * gridDim.x)

/* 获取线程ID：1D块和2D网格 */
// 计算全局线程索引，公式：全局线程ID = 块内线程索引 + 块ID * 每块线程数
#define get_tid() (blockDim.x * get_bid() + threadIdx.x)

/* 预热GPU函数声明 */
void warmup();
/* 获取时间戳函数声明 */
double get_time(void);
/* 主机端向量加法函数声明 */
void vec_add_host(FLOAT *x, FLOAT *y, FLOAT *z, int N);

//---------------------------------------------------------------------
/* 设备端向量加法内核函数（在GPU上执行） */
__global__ void vec_add(FLOAT *x, FLOAT *y, FLOAT *z, int N) {
    int idx = get_tid();  // 计算当前线程的全局索引
    if (idx < N)          // 边界检查，防止越界
        z[idx] = x[idx] + y[idx] + z[idx]; // 执行计算：z = x + y + z（注意此处非常规加法）
}

//---------------------------------------------------------------------
/* 主机端向量加法实现（在CPU上执行） */
void vec_add_host(FLOAT *x, FLOAT *y, FLOAT *z, int N) {
    for (int i = 0; i < N; i++) {
        z[i] = x[i] + y[i] + z[i]; // 与GPU内核相同的计算逻辑
    }
}

//---------------------------------------------------------------------
/* 时间获取函数实现 */
#if USE_UNIX // Unix/Linux/MacOS时间获取方案
#include <sys/time.h> // 包含gettimeofday函数头文件
double get_time(void) {
    struct timeval tv;                 // 时间结构体（秒+微秒）
    gettimeofday(&tv, NULL);           // 获取当前精确时间
    return tv.tv_sec + (double)tv.tv_usec * 1e-6; // 转换为秒为单位
}
#else // Windows时间获取方案
#include <windows.h>
double get_time(void) {
    LARGE_INTEGER timer, freq;         // 高精度计时器结构体
    static int init = 0;
    if (!init) QueryPerformanceFrequency(&freq), init = 1; // 初始化频率
    QueryPerformanceCounter(&timer);   // 获取当前计时器值
    return (double)timer.QuadPart / freq.QuadPart; // 转换为秒
}
#endif

//---------------------------------------------------------------------
/* 预热内核（空操作，仅用于初始化GPU上下文） */
__global__ void warmup_knl() {
    int i = 1, j = 2;
    i = i + j; // 无实际作用的计算，仅为触发GPU初始化
}

/* 预热GPU（避免首次启动内核的冷启动延迟） */
void warmup() {
    for (int i = 0; i < 8; i++) {      // 多次启动内核确保初始化
        warmup_knl<<<1, 256>>>();      // 启动1个块，每块256线程
        cudaDeviceSynchronize();      // 同步等待内核完成
    }
}

//---------------------------------------------------------------------
int main() {
    int N = 20000000;                   // 数据总量：2千万个元素
    size_t nbytes = N * sizeof(FLOAT);  // 计算总字节数

    /* GPU内存指针 */
    FLOAT *dx = NULL, *dy = NULL, *dz = NULL;

    /* CPU内存指针 */
    FLOAT *hx = NULL, *hy = NULL, *hz = NULL;

    /* 块大小设置（1D结构，每块256线程） */
    int bs = 256; 
    /* 计算2D网格维度 */
    // 计算需要的块数：(N + bs -1)/bs 向上取整，再开平方得到二维网格边长
    int s = (int)ceil(sqrt((N + bs - 1) / (double)bs));
    dim3 grid(s, s); // 定义二维网格布局（s x s个块）

    int itr = 30;     // 迭代次数（减少计时误差）
    double th, td;    // 时间变量（host_time, device_time）

    warmup(); // 预热GPU，避免冷启动误差

//     cudaMalloc((void**)&dx, nbytes); 是 CUDA 中用于分配 GPU 显存的关键函数调用，这里涉及到双指针（二级指针）的概念。我会用通俗易懂的方式逐步解释：
// 1. 先说基础原理
// 假设我们有一个指针变量 FLOAT* dx（指向 float 类型数据的指针），现在要在 GPU 上分配内存，并让 dx 指向这块新分配的内存。
// 但是 C/C++ 语言中函数的参数传递是值传递的，如果直接传递指针 dx，函数内部只能修改指针指向的内容，无法修改指针本身的地址值。
//因此需要用指针的指针（二级指针）来让函数可以修改指针变量的值。

// 2. 函数原型分析
// cudaMalloc 的函数原型如下：
// cudaError_t cudaMalloc(void** devPtr, size_t size);
// 参数 devPtr 是一个二级指针（void** 类型）
// 参数 size 是要分配的内存大小
// 函数功能：在 GPU 上分配 size 字节的内存，并将分配的内存地址写入 devPtr 指向的指针变量

// 双指针的作用：让函数可以修改指针变量的值（即让指针指向新分配的内存）
// 类型转换 (void**)：为了匹配 CUDA 函数原型的要求
// 最终效果：dx 指向 GPU 上的内存，后续可以通过 dx 操作 GPU 数据

    /* 分配GPU内存 */
    cudaMalloc((void**)&dx, nbytes);
    cudaMalloc((void**)&dy, nbytes);
    cudaMalloc((void**)&dz, nbytes);
    
    if (!dx || !dy || !dz) { // 检查分配是否成功
        printf("GPU内存分配失败\n");
        return -1;
    }
    printf("GPU分配: %.2f MB\n", nbytes/(1024.0f*1024.0f)); // 打印显存占用

    /* 分配CPU内存 */
    hx = (FLOAT*)malloc(nbytes);
    hy = (FLOAT*)malloc(nbytes);
    hz = (FLOAT*)malloc(nbytes);
    if (!hx || !hy || !hz) {
        printf("CPU内存分配失败\n");
        return -2;
    }
    printf("CPU分配: %.2f MB\n", nbytes/(1024.0f*1024.0f)); // 打印内存占用

    /* 初始化数据（全部赋值为1） */
    for (int i = 0; i < N; i++) {
        hx[i] = 1;
        hy[i] = 1;
        hz[i] = 1;
    }

    /* 拷贝数据到GPU */
    cudaMemcpy(dx, hx, nbytes, cudaMemcpyHostToDevice);
    cudaMemcpy(dy, hy, nbytes, cudaMemcpyHostToDevice);
    cudaMemcpy(dz, hz, nbytes, cudaMemcpyHostToDevice);

    /* GPU计算（迭代30次） */
    cudaDeviceSynchronize(); // 确保之前操作完成
    td = get_time(); // 开始计时
    for (int i = 0; i < itr; i++) {
        vec_add<<<grid, bs>>>(dx, dy, dz, N); // 启动内核
    }
    cudaDeviceSynchronize(); // 等待所有内核完成
    td = get_time() - td;    // 计算总耗时

    /* CPU计算（迭代30次） */
    th = get_time(); // 开始计时
    for (int i = 0; i < itr; i++) {
        vec_add_host(hx, hy, hz, N); // 执行主机端加法
    }
    th = get_time() - th;    // 计算总耗时

    /* 输出结果（保留科学计数法，显示加速比） */
    printf("GPU时间: %.3e秒, CPU时间: %.3e秒, 加速比: %.2f倍\n", td, th, th / td);

    /* 释放资源 */
    cudaFree(dx); // 释放GPU显存
    cudaFree(dy);
    cudaFree(dz);
    free(hx);     // 释放CPU内存
    free(hy);
    free(hz);

    return 0;
}